library(lars)
library(energy)
library(DAAG)
X=data(diabetes)
attach(diabetes)
diabetes2<-diabetes
rm(diabetes)

#begin analysis here

#thihs is for primary and interactions
#x=diabetes2$x2

#primary only
x=diabetes2$x
y=diabetes2$y

#fit the full model (include everything)
full.model = lm(y~x)
full.model.residuals = full.model$residuals
#test to see if interactions would be necessary
interactions.necessary.test = dcov.test(full.model.residuals,x,R=1999)
interactions.necessary.pval = interactions.necessary.test$p.value
interactions.necessary.pvals = c(interactions.necessary.pvals,interactions.necessary.pval

#start lars
lars.object=lars(x, y, type = c("lar"),trace = FALSE, normalize = TRUE, intercept = TRUE)
plot(lars.object)
dev.new()
plot.lars(lars.object,xvar=c("step"),plottype="Cp")

k=as.numeric(names(which.min(lars.object$Cp)))+1

summary(lars.object)
names(lars.object)
lars.orders = names(unlist(lars.object$actions))
formula.string = "y ~ "

for(i in 1:(k-1)){
	ordersplus = paste(lars.orders[i]," + ", sep = "")
	formula.string = as.character(paste(formula.string,ordersplus,sep=""))
}
formula.string = paste(formula.string,lars.orders[k],sep="")
regression.formula = formula(formula.string)
#iterate here
p.val.aggregator = NULL
lars.object$beta

x.new = as.matrix(x[,seq(1,k)])
chosen.regressors = x[,seq(1,k)]
dataframe = as.data.frame(cbind(y,x))


#ols.regression= lm(y~x.new, data = dataframe)
#ols.regression= glm(regression.formula, data = dataframe )
ols.regression1= lm(y ~ bmi + ltg + map + hdl + bmi:map + age:sex + glu^2 + bmi^2 + age:map + age:glu + sex + glu + age:ltg + age^2 + sex:map + map:hdl, data = dataframe)
summary(ols.regression1)
ols.regression2= lm(y ~ bmi + ltg + map + hdl + bmi:map + age:sex + glu^2 + bmi^2 + age:map + age:glu + sex + glu, data = dataframe )
summary(ols.regression2)

y.covtest.var = length(ols.regression$residuals
dev.new()
plot(ols.regression1$fitted.values,ols.regression2$fitted.values)
x.covtest.var = x[,seq(k+1,length(names(x)))]
#for(i in 1:100){

#test for independence
dcov.test = dcov.test(y.covtest.var,x.covtest.var,R=1999)
dcov.test.pvalue=dcov.test$p.value
dcov.test.pvalue
p.val.aggregator=c(p.val.aggregator,dcov.test.pvalue)
#}
#try 1-less k...
#k=k-1


names(ols.regression)
names(dataframe)
new.names = c("y",lars.orders[1:k])
new.indices=match(new.names,names(dataframe))
new.dataframe=dataframe[,new.indices]
names(new.dataframe)
names(new.dataframe)
cross.val.object1 = cv.glm(dataframe , ols.regression, K=10)
cross.val.object2 = cv.lm(df=dataframe ,seed=5000, form.lm=formula(y ~ bmi + ltg + map + hdl + bmi:map + age:sex + glu^2 + bmi^2 + age:map + age:glu + sex + glu), m=2,printit=TRUE)
cross.val.object3 = cv.lm(df=dataframe ,seed=5000, form.lm=formula(y ~ bmi), m=2,printit=TRUE)
names(cross.val.object)
cross.val.object$ss
k.16.ss= cross.val.object1[1]
k.16.df= cross.val.object1[2]

k.16.ss/k.16.df

cross.val.object

names(new.dataframe)

names(cross.val.object1)

#run this last...
detach(diabetes2)
